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ABSTRACT 

The ambient hot intra-halo gas in clusters of galaxies is constantly fed and stirred by 
in-falling galaxies, a process that can be studied in detail with cosmological hydro- 
dynamical simulations. However, different numerical methods yield discrepant predic- 
tions for crucial hydrodynamical processes, leading for example to different entropy 
profiles in clusters of galaxies. In particular, the widely used Lagrangian smoothed 
particle hydrodynamics (SPH) scheme is suspected to strongly damp fluid instabili- 
ties and turbulence, which are both crucial to establish the thermodynamic structure 
of clusters. In this study, we test to which extent our recently developed Voronoi 
particle hydrodynamics (VPH) scheme yields different results for the stripping of gas 
out of infalling galaxies, and for the bulk gas properties of cluster. We consider both 
the evolution of isolated galaxy models that are exposed to a stream of intra cluster 
medium or are dropped into cluster models, as well as non-radiative cosmological sim- 
ulations of cluster formation. We also compare our particle-based method with results 
obtained with a fundamentally different discretisation approach as implemented in 
the moving-mesh code AREPO. We find that VPH leads to noticeably faster stripping 
of gas out of galaxies than SPH, in better agreement with the mesh-code than with 
SPH. We show that despite the fact that VPH in its present form is not as accurate 
as the moving mesh code in our investigated cases, its improved accuracy of gradient 
estimates makes VPH an attractive alternative to SPH. 
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1 INTRODUCTION 

In the hierarchic al structure formation model 
(jWhite Frenklll99lh . interactions or mergers of galaxies 
with other galaxies are common processes, and depending 
on the environment, a diverse set of physical processes can 
become important. For example, when a galaxy falls into a 
cluster of galaxies, gravity may affect its integrity through 
gravitational tidal forces. The closer the galaxy gets to 
the cluster the more it is exposed to the hot intra-cluster 
medium (ICM). This is experienced as a headwind by the 
galaxy that compres ses and removes the gas of the galaxy 
Kjrunn GottlfioTi ). This may lead to the formation of a 
bow shock in front of the galaxy if it supersonically ploughs 
through the ICM. 

Whereas the gas at the front of the galaxy is pri- 
marily compressed by ram pressure in this situation, the 
sides are exposed to strong shear flows which may trig- 
ger fluid instabilities, leading to stripping and mixing of 



the gas with the ICM. This in turn crucially determines 
how rapidly star formation is turned off in an infalling 
galaxy, and how m etals are mixed into the group or clus- 
ter ICM (see e g. iLarson et al.1 llQSrt 'Balog h et all I2OOOI : 
IPomainko et al.l I2OO6I ). The hydrodynamical interaction 
processes between galaxies and groups/clusters may thus in- 
fluence important ob servational effects, such as the density 
morphology relation (iDresslerlllQsA Ivan der Wei et al"]l201Ql : 
iBoselli Gavazzill2006fL or the existence of a tight red clus- 
ter g alaxy sequence (|Baldrv et al ] |2004l : ICortese Hughei 
I2OO9I I. Another interesting aspect besides the stripping is 
the generation of turbulence by galaxies infalling into clus- 
ters. Numerical simulations ( Dolag et al. I2OO5I: IVazza et al] 



I2OO6I : llapichino &; Niemeved [20081 : iDolag et al .' '2009^ indi- 
cate that turbulence can be generated in the wake of an in- 
falling galaxy. This interesting phenomenon may play an im- 
portant role in determining the dynamics of the cluster gas 
and could in principle be detected through X-ray observa- 
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tions of the broadening of sharp metal hnes (|Sunvaev et al.l 
l2QQ3l ). 

Modelhng thes e processes rehably goes back to 
iGunn GottI (|l972l ) and is a considerable challenge, as 
direct hydrodynamic simulations of the relevant processes 
are very difficult to model due to the large dynamic range, 
and the uncertainties associated with the modelling of star 
formation and feedback processes. Furthermore, it is not 
clear whether the hydrodynamical techniques presently in 
use are strongly affecte d by systematic errors in this regime 
(e.g. Agertz et al. '2007). Indeed, previous studies of galaxy- 
cluster interactions with different numerical schemes dis- 
agreed strongly about how fast a galaxy loses i ts gas to the 
cluste r. For example, a grid-based simulation of lQuilis et al.l 
(|2000l ) predict ed a m uch higher stripping rate than SPH- 
based work bv Ubadi et al.1 ((1999). 

In this paper, we concentrate on the numerical as- 
pects of the stripping of a galaxy's gas as it interacts 
with the ICM. "Agertz et ajj (jioOT. ) has recently shown that 
the popular smoothed particle hydrodynamics (SPH) tech- 
nique has problems to properly account for fluid instabil- 
ities, prompting significant concerns about a possibl e un- 
phys i cal suppression of stripping processes (see also [Prio 



l2008VWadslev et al."2008; Read et al. 2010: ' Cha et~al 
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,Valcke et al. 2010; Junk et al. 2010; Abel 20^. In a recent 
study (|Hefi Springell [20101) we have therefore proposed 
a new 'Voronoi particle hydrodynamics' (VPH) m ethod 
that improves on the widely used SPH tec hnique (|LucvI 
Il977l : iGingold Monaghanl l 19771 : lLarsonlll97 8^ in several re- 
spects. By employing a Voronoi tessellation for the local den- 
sity estimate, a consistent decomposition of the simulation 
volume is achieved, contact discontinuities can be resolved 
much more sharply, and a 'surface tension' effect across them 
is avoided. Our prelim inary tests of VPH based on the 'blob- 
test' of Agertz et al.l (^2007 ) already suggested that stripping 
is more efficient in VPH compared with SPH. Here we shall 
investigate this in more detail using more realistic set-ups 
that mimic galaxy evolution processes. In addition to the 
two particle-based Lagrangian schemes for hydrodynamics, 
SPH and VPH, we will also carry out co mparison simula - 
tions with the moving- mesh code AREPO (|SDringelll2010al ). 
Whereas AREPO uses a Voronoi tessellation as well, this 
code employs an entirely different methodology for fluid dy- 
namics, based on a finite volume Godunov scheme that cal- 
culates hydrodynamical fiuxes with a Riemann solver across 
mesh boundaries. The comparison of this diverse set of three 
numerical methods is useful to understand and quantify the 
systematic uncertainties of the different methods. 

In our tests simulations, we ideally want to investigate 
as realistic conditions as possible, including a full treatment 
of gravity and dark matter, as well as star formation and 
cooling. This in principle calls for simulations in which a 
well-resolved galaxy model is dropped into an equally well 
represented cluster of galaxies. Because the substantial com- 
putational cost of such a set-up severely limits the resolution 
that can be achieved, we base part of our study on a num- 
ber of more idealised simulations, for example by placing 
galaxy models into a 'wind tunnel' that mimics the imping- 
ing fiow of gas onto a galaxy in orbit in a cluster. Finally, 
we also carry out cosmological simulations of the formation 
of the 'Santa Barbara cluster' (Frenk et al. 1999), primarily 
to study how well VPH performs in non-radiative cosmolog- 



ical simulations of cluster formation compared to the other 
techniques. Already in past studies, the Santa Barbara clus- 
ter comparison project has led to important insights into 
how numerical effects impact the thermodynamic structure 
of simulated clusters, and for this test problem, there is al- 
ready a considerable body o f results in the literature (e.g. 
Wads lev et gH l2004l : ISuring el 2005; Thacker CouchmanI 
|2006|). 

This paper is structured as follows. In Section [21 we in- 
troduce the numerical methods we will use in our numerical 
comparison study. In Section [3l we check how well the differ- 
ent numerical schemes represent the evolution of an isolated 
galaxy, and whether there are already significant differences 
at this level. We then consider in Section [4] wind-tunnel ex- 
periments in which we expose the isolated galaxy models 
to a supersonic wind, allowing a detailed examination of the 
stripping process. In Section [5l we follow up on these experi- 
ments by studying the behaviour of a galaxy model that falls 
into an isolated spherical cluster, comparing the VPH results 
with those obtained with SPH and AREPO. Finally, in Sec- 
tion [6] we investigate the performance of VPH in cosmolog- 
ical simulations of cluster formation, using re-simulations of 
rich clusters identified in the Millennium Simulation as well 
as the Santa Barbara cluster initial conditions. We give a dis- 
cussion of our results in the context of a vortex test problem 
in Section [3 combined with a summary of our conclusions. 



2 METHODOLOGY 

We begin by introducing the different codes and compu- 
tational methods we use in this study. In the interests of 
brevity, we shall only describe the most important charac- 
teristics o£_the codes; full details can be found els ewhere 
(|SDringel[[2005l : [HeB SDrin^[2010l : [SpringejlioToil ). 



2.1 Smoothed particle hydrodynamics 

Smoothed particle hydrodynamics (SPH) is a particle based 
approach to fiuid dynamics that has seen widespread use 
in astronomy since it was first introduced more than three 
decades ago (|Lucv[[1977[ : [Gingold Monaghan[[l977[ l. SPH 
discretises the mass of the fiuid in terms of particles, whose 
dynamics is governed by a fiuid Lagrangian of the form 



^ ' ^rriivl - rrii Ui{pi, Si) 



(1) 



Here Ui describes the thermal energy per unit mass, which 
for an ideal gas depends only on density pi and specific 
entropy Si. An adaptive, spherically symmetric smoothing 
kernel is employed to estimate the density based on the spa- 
tial distribution of an approximately fixed number of near- 
est neighbours. (As a standard value within this work we 
set the number of neighbours to N — 48.) The Lagrangian 
in Eq. ([!]) then uniquely determines the equations of mo- 
tion in a form that simultaneously conserves e nergy and en- 
tropy (Springel k Her nauist'' 2002l : [Price k Monaghanl l2007l : 
.Rosswog 2009 : Springel 20101J ). However, an artificial viscos- 
ity needs to be added to allow for dissipative shock capturing 
and to produce a damping of post-shock oscillations. If not 
explicitly denoted otherwise we use the viscosity parametri- 
sation of the GADGET code with a viscosity strength of 
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a = 1, and we apply the viscosity limiter in the presence of 
shear introduced by Balsara (1995). 

Among the primary advantages of SPH are its very good 
conservation properties, the manifest Gahlean invariance at 
the discretised level of the equations, and the absence of 
preferred spatial directions. However, it has recently become 
clear that contact discontinuities, in particular, are challeng- 
ing for the method. For example, they are associated w ith a 
spurious surface tension effect in SPH (|SDringelll2010bh . and 
fluid instabilities such as t he Kelvin- Helmholt z instability 
are suppressed across them (|Agertz et al.ll200'7l V 

2.2 Voronoi particle hydrodynamics 

Because of the accuracy problems of SPH in certain situ- 
ations, we have recently developed another particle-based 
hydrodynamic method that differs from SPH in impor- 
tant ways, while still sharing a number of similarities. 
In this ^Vor onoi Particle Hydrodynamics' (VPH) method 
(|HeB Spr ingel 2010), the fluid is also discretised in terms 
of mass elements m^, and the same fluid Lagrangian as in 
Eq. ([1]) is used. Also identically to SPH we assume a poly- 
tropic equation of state P — Ap^ with an adiabatic co- 
efficient 7 in all our tests. However, the kernel estimation 
technique of SPH is replaced by a very different approach 
to calculate the gas density. This is done by constructing a 
Voronoi tessellation for the point set, in which to each of the 
particles the volume of all the space which is closer to this 
point than to any other point is assigned. The density can 
then be obtained as the particle mass divided by the volume 
of the associated Voronoi cell. We note that such a tessel- 
lation technique allows a close to optimum exploitation of 
the d ensity information con tained in the particle distribu- 
tion (Pelup essv et all 12003 !). and in particular, very sharp 
discontinuities can be resolved over the distance of a single 
ceh. 

With the density computed, and given specified en- 
tropies Si for every particle, we can use the ideal gas La- 
grangian to derive unique equations of motion for VPH. 
The key quantity that is needed in doing this is the deriva- 
tive of the particle volume with respect to the coordinates 
of the surrounding particles, but thanks to the mathemati- 
cal properties of the Voronoi tessellation, the correspond- 
i ng geometric factors can be computed relati vely easily 
(ISerrano Espanoll lioOll: I kefi Springej linol) and have 
a clear geometric meaning. For example, the primary force 
between two VPH particles is their pressure difference times 
the area of their common tessellation face. Because the dy- 
namics in VPH is derived from the Lagrangian given in 
Eqn. energy, momentum and entropy are conserved ex- 
actly, just like in SPH. However, there is no surface ten- 
sion effect, and the VPH particle volumina add up to the 
total simulated volume exactly, something that is not the 
case in SPH in general. One can hence hope that VPH 
provides an inte resting improvement over SPH. Indeed, in 
iHeB k Sprint (2010) we have shown that this is the case 
for several test problems, such as the growth of Kelvin- 
Helmholtz instabilities at contact discontinuities with a large 
density jump. It is one of the primary goals of this paper to 
see whether these improvements are also relevant in cosmo- 
logical structure formation problems. 

In VPH, we need to rely on an artificial viscosity to 



treat shocks, simi larly to SPH. We here follow the standard 
SPH approach bv lGingold MonaghanI (p-977) and Bal saral 
(1995) in introducing an extra friction force that reduces the 
kinetic energy and transforms it into heat in regions of rapid 
compression. Due to its similarities to its SPH counterpart 
we usually use a value of a = 1 for our artificial viscosity in 
VPH as well, but we also vary a explicitly to demonstrate 
its implications. In addition, in VPH we can also add some 
weak, non-dissipative "shape forces" that are designed to 
maintain a reasonably regular mesh geometry, which tends 
to improve the overall accuracy of the VPH d ynamics. As 
discussed extensively in HeB & S^rm^ ([joTo), there are 
different possibilities to introduce such extra forces. We de- 
rive them by adding extra terms to the Lagrangian which 
penalise highly distorted mesh cells. For the shape correction 
scheme we usually choose /3o = 0.2 and /3i = 0.01. However, 
the optimum choice for the strength of these terms is unfor- 
tunately not (yet) clear. For this reason, we have typically 
carried out all our tests with several different settings for the 
strength of the artificial viscosity and the shape forces, find- 
ing in most cases that our results are quite insensitive to the 
detailed choices over a broad range of parameter settings. 

Both the SPH and VPH simulations we study here have 
have been carried out with GADGET-3, an improved version 
of the publicly availa ble co smological code GADGET-2 (last 
described in SDrinid"'2005). The code is parallelized for dis- 
tributed memory machines and uses a hierarchical tree algo- 
rithm for calculation the self- gravity of the gas. A collision- 
less stellar or dark matter component can be optionally in- 
cluded as well. 

2.3 Hydrodynamical moving-mesh code 

Finally, the third numerical technique for hydrodynamics we 
examine i s impl emented in the moving-mesh code AREPO 
(|Springell l2010al V It uses an entirely different approach in 
which the volume is discretised. The hydrodynamical equa- 
tions are solved in terms of a Godunov scheme where at 
each mesh interface a Riemann problem is calculated and 
the resulting fluxes are used for an exchange of conserved 
fluid quantities between the cells. For higher spatial accu- 
racy, this is combined with a piece- wise linear reconstruction 
of the gas distribution and a second-order time integration 
scheme, yielding overall second-order integration accuracy 
in smooth parts of the flow. 

The particular mesh used by AREPO is given by the 
Voronoi tessellation of a set of mesh- generating points. The 
use of this particular type of mesh allows the method to be 
formulated such that the mesh can move along with the gas, 
continuously adjusting its topology to the fluid motion. The 
method exploits the property of a Voronoi tessellation that 
the mesh changes continuously when the generating points 
are moved, without showing problematic mesh-tangling ef- 
fects. If the points are moved with the local flow velocity 
(which is the standard way to operate the code) , the method 
then has a similar Lagrangian character as SPH/VPH. How- 
ever, unlike in SPH and VPH, the masses of AREPO's resolu- 
tion elements (the cells) are not required to remain strictly 
constant. In fact, if the mesh- generating points are being 
fixed on a static Eulerian grid, then AREPO is equivalent to 
the MUSCL-Hancock scheme, a popular method for second- 
order Eulerian gas dynamics on structured grids. AREPO 
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should hence be viewed as being closely related to ordinary 
Eulerian mesh-based hydrodynamics, except that its ability 
to move the mesh along with the flow considerably reduces 
advection errors. 



3 ISOLATED GALAXIES AND THEIR 
EVOLUTION 

Before we analyse the interaction of individual galaxies 
within a cluster environment, we here study galaxy models 
in isolation to see how well the different numerical schemes 
treat them and whether our initial conditions are sufficiently 
stable. This is also intended to identify possible subtle sys- 
tematic differences between the numerical schemes that may 
be difficult to detect in more complicated set-ups. We hence 
use identical initial conditions for the three codes such that 
differences introduced by varying start-up procedure are 
minimised. We we will see however that even for identical 
particle positions, masses, and temperatures, there can be 
variations in density (and hence specific entropy since we 
start with a prescribed temperature) directly after the start 
of a simulation, as a result of the different density estimation 
methods. We are in fact especially interested in the question 
whether these small differences have any bearing on the fur- 
ther evolution of isolated galaxy models. 

To set-up compound galaxies in isolation, we con- 
struct dark matter halo models and stellar disks in dy- 
namical equilibriu m through an approximate solution of the 
Jeans equations ([HernauistI [l993; Springel & Whit^ Il999l : 
ISpringel et al.ll2005l l. The primary assumption made is that 
the local velocity dispersion can be approximated reason- 
ably well with a triaxial Gaussian. For the models considered 
here, a dark matter halo with a Hernquist profile is used, but 
alternatively a truncated NFW profile could be used as well, 
with no difference to our conclusions. The baryonic matter 
is represented with a central stellar bulge, an exponential 
stellar disk, and an exponential gaseous disk. In addition, 
we adopt a gas distribution of very low density in the halo, 
taken to be in hydrostatic equilibrium in the total potential 
such that this gas is close to the virial temperature. This 
component helps to avoid problems in VPH and AREPO 
due to "empty space" - unlike SPH, these codes can not 
easily deal with vacuum boundary conditions as they tessel- 
late space. By having the low density background gas in the 
halo, which we extend to a box region around the galaxy, 
we can avoid this problem. 

F ollowi ng the formalism and nomenclature oflMo etaL 
(|l998l ) and ISpringel et al.1 (|2005h . we have chosen a virial 
velocity of V200 = 207kms~^, yielding a total mass of 
Mtotai = V"2oo/(10 GHo) = 2.1 X 10^^ h'^MQ for the system, 
of which a fraction Mdisk /Mtotai = 0.041 is assumed to be in 
the disk, a fraction Mbuige/Mtotai = 0.01367 of the mass is 
in the bulge and the rest is distributed in a halo with a con- 
centration of c = 9.0 and a spin parameter of A = 0.033. The 
Hubble constant is set to h = ifo/(100 kms"^Mpc~^) = 0.7. 
Since we want to concentrate on the gas dynamics, we 
make the disk very gas-rich by giving it a gas fraction of 
50%. The gas component in the halo is given a mass frac- 
tion Mgas, halo /Mtotai = 0.0025, Considerably less than the 
amount of gas in the disk. This additional gas in the halo is 
represented with particles/cells of mass equal to the gas par- 
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Figure 1. Vertical gas density profile of our isolated galaxy model 
at the initial time t = when simulated with different numerical 
techniques. The green dashed lines gives the actual initial density 
profile obtained by averaging the mass distribution over the radial 
range 1.5/i~^kpc < R < 12/i~^kpc. The black lines correspond 
to the mean density estimated for the particles falling in the cor- 
responding region in our SPH simulation. The blue line is the 
corresponding measurement for VPH. Note that the actual par- 
ticle distribution is identical for SPH and VPH at time t = 0, so 
the difference merely refiects different systematics in the density 
estimation techniques employed by the two codes. We also note 
that the result for AREPO is the same as the one for VPH, since 
they both use the same initial conditions and base their density 
estimate on the Voronoi tessellation. 



tides in the disk, and its spatial distribution is clipped out- 
side a box of side- length 2000 h~^kpc centred on the galaxy, 
allowing us to impose periodic boundary conditions for the 
gas there. 

We include radiative cooling (for a primordial com- 
position), and model star formation and supernova feed- 
back with the simp le multi-phase model for the ISM of 
ISpringel HernauistI ([20031). For the supernova feedback 
model we adopt an elevated effective temperature of Tsn = 
lO^K and an evaporation efficiency of Aq = 10^. This is 
done in order to shift the density threshold for the onset of 
star-formation to a 10 times higher v alues than in the de- 
fault model of ISpringel HernauistI (|2003h , while keeping 
the star- format ion timescale at — 2.1 Gyr and maintain- 
ing a good agreement with the iKennicuttI (|l998h relation. 
With this change, the star formation threshold corresponds 
to a hydrogen number density of nu = 1.3 cm~^. Such an 
elevated star formation density threshold has recently been 
advocate d to facilitate the fo rmation of more realistic spiral 
galaxies (|Guedes et al.|[201ll V and it will sharpen any po- 
tential issues brought about by the density contrast between 
disk and gaseous corona, which is what we are specifically in- 
terested in in this work. If not denoted otherwise we employ 
this model throughout this paper whenever star formation 
and/or cooling is required. 

In our default intermediate resolution (hereafter called 
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Figure 2. Vertical gas density profiles and sampling point densi- 
ties at an evolved time {t = 0.25 /i~-'^Gyr) for an isolated galaxy 
model, simulated either with SPH (black), VPH (blue) or AREPO 
(red). The top panel gives the gas density profile obtained by av- 
eraging the density estimates obtained for each simulation parti- 
cle/cell as solid lines, while the dashed lines are the actual mean 
mass profile by summing up the masses and dividing by the bin 
volumes. We see that for VPH and AREPO these measures agree 
reasonably well with each other, while in SPH some larger sys- 
tematics at the strong contact discontinuity between cold/dense 
disk and hot /tenuous corona appear. In particular, there is in 
fact a small region with essentially no tracer particles in SPH, 
which is more explicitly seen in the bottom panel, which sim- 
ply gives the total number of tracer particles (fiuid particles or 
mesh-generating points, respectively) found in the corresponding 
region. 
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'R4'), the gas in the disk is represented by 80000 particles 
and an equal number of collision-less particles for the stellar 
disk. The bulge is represented by 40000 particles, and the 
dark matter halo by 120000 heavier particles that are given 
a larger gravitational softening length of 0.5/i~^kpc com- 
pared to 0.25 /i~^kpc for the softening length of the baryonic 
particles. The gas in the halo out to the boundaries of the 
simulation box is represented with an additional 9756 parti- 
cles. Besides this default resolution we have also simulated 
realisations of the same galaxy model both at lower and 
higher resolutions, in order to get a good sense of numerical 
resolution effects. To this end we lowered the resolution in 
three steps by factors of 2 in all particle components, and 
we increased it in three steps by factors of 2 as well. We 
hence arrived at 7 different resolutions R1-R7, with 10000, 
20000, 40000, 80000, 160000, 320000, and 640000 resolution 
elements in the gaseous disk, respectively (see Table [1]). The 
gravitational softening was varied in proportion to the cube 
root of the mass resolution, e oc m^^^, in these simulations. 
All our simulations with SPH, VPH and the AREPO code 
were started from the same identical initial conditions files in 
order to yield a direct comparison that also allows an identi- 
fication of different start-up systematics. For the simulations 
with AREPO, we have a ctivated on-the-fiy refinem ents and 
derefinements similar to IVogelsberger et al] (|201ll ) in order 
to keep the mass resolution always close to the initial value. 

3.1 Gas density maps and structure of the disk 

As we have remarked above, identical particle positions and 
masses do not necessarily imply equal densities in our dif- 
ferent schemes. To examine this point, we consider in Fig- 
ure [1] each method's estimate of the initial gas density pro- 
file in the ^-direction, perpendicular to the disc. The solid 
curves represent averages of the densities estimated for the 
individual particles/cells as a function of distance from the 
disk plane, while the dashed line is the actual average mass 
density distribution obtained by binning the particles di- 
rectly. Since both VPH and AREPO use a Voronoi tessella- 
tion for the density estimate, they agree on the densities at 
the initial time, but the kernel estimates of SPH (based on 
^ngb = 48 smoothing neighbours) show an interesting dif- 
ference. In particular, the particles close to the dense disk 
gas (for z ^ 0.2 — 2.0 /i~^kpc) show an elevated density esti- 
mate because they "see" some of the dense particles within 
their smoothing radius. This effect is not unexpected given 
that sharp discontinuities will always be smoothed out to 
some extent in SPH, by construction. The density estimate 
of VPH has a somewhat better resolving power, but here a 
different systematic effect becomes noticeable. Particles in 
the surface layer of the dense disk sometimes extend their 
volume quite far into the region above the disk until they 
"see" a particle of the background corona, thereby under- 
estimating the density in this region. Finally, there is an 
interesting systematic difference between SPH and VPH at 
the very centre of the disk, where SPH shows a small posi- 
tive bias in the density estimate. This comes about in part 
due to the Poisson sampling present in the initial conditions, 
which in SPH causes larger positive excursions in the density 
estimates. For a system in pressure equilibrium, this effect 
is weaker, but it not necessarily vanishes because also here 
the SPH density may be biased and the sum of the volumes 
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run name i 


gas particles 


disk -|- bulge particles 


softening gas/stars 


DM particles 


softening DM 


rti 


iUUUU 


1 nnnn i c^nnn 
iUUUU + oUUU 


U.o ri Rpc 


1 Knnn 


i.U Al Kpc 


R2 


20000 


20000 + 10000 


0.4 h-^kpc 


30000 


0.8 h-^kpc 


R3 


40000 


40000 + 20000 


0.31 /i-^kpc 


60000 


0.6 h-^kpc 


R4 


80000 


80000 + 40000 


0.25 h-^kpc 


120000 


0.5 h-^kpc 


R5 


160000 


160000 + 80000 


0.2 /i-^kpc 


240000 


0.4 h-^kpc 


R6 


320000 


320000 + 160000 


0.16 /i-^kpc 


480000 


0.32 h-^kpc 


R7 


640000 


640000 + 320000 


0.125 h-^kpc 


960000 


0.25 h-^kpc 



Table 1. List of numerical parameters for the simulations shown in the resolution study of Figure [6] The columns give the simulation 
name, the initial number of gas particles, the initial number of star particles in disk and bulge, the gravitational softening lengths for the 
gas particles, the number of dark matter particles, and finally the gravitational softening length of the collisionless dark matter particles 
in the halo. 



associated with each particle is not guaranteed to add up to 
the total volume. 

In Figure [2l we analyse the vertical structure of the 
gas after the isolated galaxy has been evolved for a time 
of 0.25/i~^Gyr. In the top panel of the figure, we again 
compare measurements of the mean density estimated for 
particles/cells as a function a distance z above the disk 
mid-plane with the actual mean mass profile present in the 
simulation. We give results for all three numerical methods 
(solid lines), and since the actual gas distribution can have 
evolved differently in the three methods, this is given indi- 
vidually as well (dashed lines). We see that the two measures 
of the density stratification track each other well both in 
VPH and AREPO. However, the moving- mesh code shows 
higher density just above the disk compared with VPH. If 
anything, SPH shows still higher densities in this region, al- 
though its actual mass distribution is very similar to that 
of VPH. This situation arises due to a significant difference 
SPH exhibits in the estimated density versus the real mass 
distribution in the regime of the 'shoulder' of the disk's den- 
sity distribution. Also, there is actually a 'gap' in the run of 
SPH's density estimates, arising simply because we find no 
SPH tracer particles for measuring the density there. This 
gap becomes more explicit in the bottom panel of Figure (2] 
where we show counts of fluid tracer particles in logarith- 
mic bins in the same region above the disk. The sampling 
gap ai z ^ 1 h~^kpc in SPH is evident in this figure , but it 
does not occur in VPH. We note that lAgertz et all (l2007l ) 
pointed out that such a gap is seen generically at strong 
density jumps in SPH, and that it may be a primary cause 
for the suppression of fluid instabilities across s uch contact 
discontinuiti es. Recently, a number of authors ( Pried [ioosi : 
iRead Hav field 2011) have suggested that SPH should be 
outfitted with additional dissipation schemes that smooth 
out such structures to improve on this behaviour. 

In Figure O we show face-on maps of the galaxy's 
projected gas density distribution, after a period of t = 
0.5/i~^Gyr of evolution. It is reassuring but probably not 
too surprising that the overall morphology is extremely sim- 
ilar, as the gravity from the dark matter halo and the stel- 
lar disk are primary driving forces of the disc dynamics. 
However, upon close inspection, one can identify some in- 
teresting systematic differences between the techniques that 
manifest themselves in the gas density maps. Compared with 
SPH, the higher effective spatial resolution of VPH produces 
crispier but marginally noisier looking features such as spi- 
ral arms. Perhaps the most prominent difference is however 



the lower density of the gas found in VPH compared to SPH 
in some regions between spiral arms. The gas distribution of 
the mesh-based AREPO looks a bit sharper and less noisy 
than both particle-based schemes. 

Nevertheless, there is reassuring similarity of the gas 
surface density distribution between the different tech- 
niques, an impression that is corroborated by maps of the 
projected stellar mass density in the disks, which we show 
in Figure |4l Here we show only the newly formed stars in 
the simulations, so that any structure in this component di- 
rectly reflects the underlying gas dynamics. We find good 
agreement for these collisionless particles, which can also 
be interpreted as evidence that the gravitational dynamics 
of collisionless particles is followed with equal quality in all 
three codes. This verifies that differences in the evolution 
of our galaxy models are expected to arise only from the 
different treatment of the hydrodynamics, and not from dif- 
ferences in the way the collisionless dynamics of stars and 
dark matter is followed. 

Despite the systematics we identified in the density es- 
timates and in the projected gas distributions, we overall 
find that all three codes evolve the galaxy model in a sim- 
ilar and stable fashion. This is further confirmed by radial 
surface density profiles of the gas and the newly formed stars 
shown in Figure [5] All three methods are clearly able to re- 
tain the initial structure of the galaxy model to a similar 
degree, yielding only very minor differences after some time. 

3.2 Star formation rate evolution 

Despite the good agreement between the different methods 
noted above, some minor but interesting discrepancies are 
revealed when looking in detail at the time evolution of the 
star formation rate at different resolutions, as done in Fig- 
ure [6l In the top panel, we compare the results for SPH, 
VPH and AREPO at our default resolution of R4, which lies 
in the middle of our extended set of models used in our reso- 
lution study. There is in general good agreement between the 
runs at late times, where SPH and AREPO agree extremely 
well and VPH lies at most a few per cent lower. However, at 
earlier times some large differences are noticeable. First of 
all, at t = 0, the star formation rate for SPH is about ~ 20 
per cent higher than that of VPH and AREPO, a difference 
entirely caused by different density estimation methods (see 
below). But after an initial transient, when the disk set- 
tles from the approximate equilibrium realised in the initial 
conditions to a proper self-consistent equilibrium, somewhat 
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Figure 3. Projected gas density maps dX t = 0.5/i~-'^Gyr of our 
isolated galaxy model simulated with SPH (top), VPH (middle) 
and AREPO (bottom). Each map measures 30 /i~-'^kpc on a side, 
employs a logarithmic colormap, and was constructed with an 
identical adaptive binning method. 
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Figure 4. Projected mass density of newly formed stars 
(i.e. omitting stars already present in the stellar disk of the ini- 
tial conditions) in a galaxy model evolved to time f = 0.5 h~^GYv 
with three different numerical schemes, SPH, VPH, and AREPO, 
from top to bottom. The maps measure 30 /i~^kpc on a side. 
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Figure 5. Surface density profiles of gas (solid lines) and newly 
formed stars (dashed lines) at time t = 0.5/i~-'^Gyr for our iso- 
lated galaxy model when simulated with either SPH (black), VPH 
(blue) or AREPO (red). 



larger differences are present for an extended period of time. 
VPH tends to show a slightly smaller star formation rate 
than SPH, with AREPO lying somewhere in the middle. 

Some clues to the origin of these differences are pro- 
vided by the lower three panels of Figure [H which give the 
evolution of the SFR for all seven numerical resolutions we 
considered for our isolated galaxy model, and for each nu- 
merical method. We see that all three techniques show some 
residual drift in the SFR at low resolution, but they tend to 
converge reasonably well towards higher resolution. Given 
the fact that the mass resolution is varied by a factor of 64 
here, a scatter in the star formation rate of order 10 — 20 per 
cent can probably be considered satisfactory. Nevertheless, 
there are clearly interesting differences in the strength of the 
resolution-dependence of the numerical techniques. In par- 
ticular, VPH shows a somewhat stronger reduction of the 
SFR at the lowest resolution compared to SPH, whereas the 
latter is remarkably resilient to resolution changes. Perhaps 
counter-intuitively, this arises despite the relatively strong 
biases in SPH's density estimate at the interface between 
disk and corona. 

Another look at this issue is provided by Figure [71 in 
which we show how much gas mass is estimated by the codes 
to lie at a certain density value. In the top panel of Fig. [71 
this is given for the initial time. The fact that SPH system- 
atically estimates more gas to lie at a density value above 
the star formation threshold than the Voronoi-tessellation 
schemes (despite an identical point distribution in this case) 
explains the offset in the SFRs at t = 0. However, also at 
later times, as seen in the bottom panel of Fig. [71 such a 
systematic difference persists. The variations in the star for- 
mation rates calculated by the different codes appear hence 
to be primarily driven by the way the sharp edge of the 
star- forming disk is represented. These differences are how- 
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Figure 6. Comparison of the time evolution of the star formation 
rate of an isolated galaxy model, simulated with SPH (black), 
VPH (blue), and AREPO (red). The top panel compares the three 
different simulation techniques directly with each other, at the 
comparatively high resolution of 160000 gas particles/cells in the 
initial disk. The other three panels show resolution studies for 
each of the codes individually, where the resolution in the gas 
was varied between 10000 and 640000 resolution elements. The 
number of dark matter and stellar bulge/disk particles has been 
varied in proportion in each of the runs of the series. We note that 
the lower resolution runs in VPH and AREPO are somewhat more 
strongly affected by resolution than for SPH, but in general all 
the methods converge and give consistent results with each other 
at high resolution. 
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Figure 7. Distribution of gas mass as a function of density esti- 
mated by our different numerical techniques (SPH in black, VPH 
in blue, and AREPO in red) in a simulation of an isolated disk 
galaxy. In the top panel, we show the distribution of gas mass in 
the initial conditions at t = 0, where the particle/cell distribution 
is identical for all three runs. Nevertheless, SPH estimates higher 
densities for some of the gas in or close to the surface of the gas 
disk, such that more gas mass ends up above the density thresh- 
old for star formation (vertical dotted line). At the later time 
of t = 0.25 h~^GyT shown in the bottom panel, the difference 
has become considerably smaller, but there is still an important 
systematic offset between VPH and SPH. 



ever quite minor overall and can be largely ignored for our 
subsequent investigations. 



4 GALAXY IN A WIND TUNNEL 

Simulating the interaction of a galaxy with the sparse intra- 
cluster gas of a galaxy clus ter at high-resolut i on is compu- 
tatio n ally challeng!:ing ( e .g. lAbadi e t al.' 'l999; Tittle v et al.l 
200l |: ISchulz Struckl l200ll : iRoed iger Briiggen |2C 



2007 , l2008h , especiallv with the particle-based methods 
SPH and VPH, which do not easily lend themselves to 
adaptive refinement techniques. In fact, they basically re- 
quire particles o f equal or very similar mass to work well 
(jOtt &: Schnetter,,20 03). Hence, if one simply wants to let a 
galaxy model fall into a massive cluster, the latter needs to 
be represented with very high particle number (because it 
is so massive), such that one ends up spending only a tiny 
fraction of the computational effort onto the galaxy and the 
surrounding gas, where the actual interaction takes place. A 
further complication is that some of the stripped gas mass 
may be distributed over a large region across the cluster, 
corroborating the need to simultaneously account for the 
whole cluster at high resolution. 



We will nevertheless consider such simulations later on 
in Section [5] to check the validity of our approach, but here 
we first investigate an alternative, more idealised set-up, 
which makes it much easier to reach an adequate resolu- 
tion. In this approach, we effectively put a galaxy model 
into a 'wind tunnel', i.e. a rectangular box in which we let 
gas stream onto the galaxy with prescribed density, veloc- 
ity and temperature, matched to what we expect during a 
cluster passage. The much smaller volume that needs to be 
simulated around the galaxy in this situation drastically low- 
ers the computational expense, and even more importantly, 
the simplification brought about by such a controlled set-up 
makes it much easier to distinguish between different nu- 
merical effects. 



4.1 Initial conditions 

For definiteness, we put our model galaxy into a rectangu- 
lar box of side- length 80/i~^kpc. At one side of the box, 
gas is constantly injected, which we realise in the case of 
SPH and VPH by creating new particles in a suitable fash- 
ion, mimicking an infinitely extended grid of particles that 
moves into the simulation domain. On the opposite side of 
the box, we effectively implement outflow boundary condi- 
tions by removing particles once they start moving out of the 
box. Since we here focus on supersonic winds, the removal of 
particles at the outflow side does not lead to the propagation 
of perturbations upstream to the inflowing gas. In the case 
of the AREPO code, the inflow and outflow regions are re- 
alised in an equivalent, yet technically different fashion. We 
here use the on-the-fly refinement and derefinement features 
of the AREPO code to produce new mesh-generating points 
in the inflow region, and to remove them near the outflow 
side. The net result is identical to the SPH and VPH cases, 
i.e. the galaxy at the centre of the box is hit by a wind of 
particles/cells with prescribed density and velocity. This im- 
pinging wind has the same gas mass resolution as our target 
galaxy model. 

For all the boundaries of the box that form the enclosing 
sides of our wind tunnel, we adopt periodic boundary condi- 
tions for the gas as far as the hydrodynamics is concerned. 
However, periodic boundary conditions are not imposed for 
self-gravity, which we fully include (note that our galaxy is 
held together by gravity). Since the mass of the gas in the 
wind is quite small, most of it is unbound to the galaxy, 
and our restricted treatment of self-gravity to the region of 
the box is still a good approximation. We note that a num- 
ber of similar win d tunnel experimen t s have previously been 
carried out (e.g. lAgertz et al.l l2007l : llapichino Niemeverl 
2008), but without a treatment of self gravity. While our 
setup is not capable of accounting for the tidal forces that 
arise when an extended galaxy travels through the gravita- 
tional potential of a cluster, these forces can be expected 
to be sub-dominant compared to the hydro dynamical forces 
in the outer parts of a cluster, which is the region we are 
most interested in. However, gravitational forces between 
dark matter, stars and gas inside the galaxy are very im- 
portant, and our approach is the first that allows a detailed 
study of the response of the gravitationally coupled disk- 
halo system to the ram-pressure of the impinging wind. In 
particular, this allows a realistic measurement of the accel- 
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Figure 8. Projected gas density maps of a galaxy (seen edge-on) that experiences a supersonic wind from below. The different panels show 
the logarithm of the density field at different times, from T = 0.04 h~^GyT (bottom row), T = 0.4 h~^Gyr (middle row) to T = 0.7 h~^GyT 
(top row), and for different numerical techniques, AREPO (left column), VPH (middle column) and SPH (right column). All maps are 
generated with the same adaptive smoothing algorithm and have an extension of 30/i~^kpc X 45/i~^kpc. 
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t [Gyr/h] 

Figure 9. Loss of star-forming ISM gas as a function of time 
due to stripping through a supersonic wind. The gas remaining 
in the ISM of the galaxy is here measured according to the cri- 
teria of Eq. ([21). The three sohd hues show the measurements for 
our simulations with SPH, VPH and AREPO, as labelled. We 
find that the stripping is quite similar in VPH and AREPO, but 
significantly less efficient in SPH. 



eration of the whole galaxy due to the ram-pressure force 
exerted by the wind. 

In our default set-up, we have chosen a constant density 
of p = 1.90 X lO^'^gcm"^ and a velocity of ^; = 3000 km s"^ 
for the wind, letting it flow onto the galaxy in a face-on ori- 
entation. A density and velocity of this magnitude would be 
typical close to apocentre for a galaxy that has fallen into a 
large cluster. Of course, in reality the strength of the wind 
will be a function of the orbital phase and of parameters such 
as cluster size and angular momentum of the galaxy orbit. 
Changing these values should however not change our results 
at a qualitative level. For the galaxy model studied in our 
simulations, we have adopted structural properties identical 
to those described in Section [S] We place the galaxy sud- 
denly into the supersonic wind, refraining from any attempt 
to impose a more gradual start-up in which the wind would 
somehow be slowly turned on. 



4.2 Stripping 

Hydrodynamically, the galaxy represents an obstacle in the 
supersonic wind, causing a bow shock in the upstream re- 
gion ahead of the galaxy. Inside the bow shock and ahead of 
the front side of the galaxy, the compressed wind is slowed 
down and exe rts sig nificant ram-pressure onto the galaxy 
(iGunn Gottill972l ). whereas around the sides, a region of 
strong shear flow is produced as the wind streams around 
the galaxy. In this region, the formation of Kelvin-Helmholtz 
instabilities is expected, which accelerate the stripping of gas 
from the galaxy's disk and mix it with the downstream flow. 




Figure 10. Maps of the projected star formation rate density in 
a galaxy model (seen edge-on) that encounters an IGM wind in 
the upwards direction. We show results at T = 0.4 h~^Gyr for our 
wind tunnel experiments, carried out with AREPO (left), VPH 
(middle) and SPH (right). Each panel extends over 36/i~-'^kpc x 
48/i-ikpc. 



lAgertz et al has previously investigated the sim- 

pler case of a spherical, non-self-gravitating gaseous obsta- 
cle (the 'blob-test') and found substantial discrepancies be- 
tween different hydrodynamical methods for the rate of gas 
stripping and the time until eventual complete disruption 
of the gas blob. Here, we are particularly interested to see 
whether such differences also occur between SPH, VPH and 
AREPO when the more realistic disk-wind setup we are in- 
vestigating is considered. 

In Figure [SI we show projected gas density maps of the 
galaxy in the wind tunnel in an edge-on orientation, with 
time evolving from bottom to top. While the morphology 
of the stripping process looks similar at early times in all 
three codes considered, significant differences develop over 
time. Both of the particle-based methods loose small clumps 
of gas and exhibit string-like features of quite dense gas. 
This effect is particularly strong in SPH. In contrast, the 
AREPO disk stays comparatively coherent, with all of the 
stripped gas moving to lower densities quickly. Hence there 
remains no dense debris in the downstream flow. At the 
latest time shown in the figure (top row), the residual gas 
disks of AREPO and VPH are visually significantly smaller 
than in SPH, suggesting that more gas has been stripped. 

This impression is born out by a quantitative measure- 
ment of the gas mass that is still in the ISM phase of the 
galaxy model. We consider a gas particle (or cell) to be part 
of the ISM of the galaxy when it is sufficiently dense and 
cold. Specifically, we require the conditions 

p > 2000 Kwind and T < 0.5 ^wind (2) 

to be met, where pwind and Twind are the density and tem- 
perature of the wind. Furthermore, we impose an additional 
condition on the maximum allowed separation r of a particle 
from the centre-of-mass of all gas that fulfils the criterion of 
Eq. (|2]). Only for r < 8kpc the gas particle is considered to 
be part of the ISM of the primary galaxy. This condition is 
introduced in order to discard cold clumps of gas that have 
been stripped out of the galaxy but have still stayed cold 
and dense. Indeed, it turns out that this secondary crite- 
rion is quite important in our SPH simulations, where the 
galactic disk tends to shed dense clumps under the action 
of the ram pressure. These fragments can still fulfil Eqn. (|2]) 
despite being stripped, but as they separate from the main 
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Figure 11. Time evolution of the sum of the mass of star-forming 
gas and the mass of all stars formed in the whole simulation box 
since the beginning of the simulations of a galaxy exposed to an 
impinging IGM wind in a wind tunnel. We show results for dif- 
ferent simulation techniques, including an SPH-based simulation 
(black solid), VPH (blue solid), and AREPO (red). 




Figure 12. Stellar density maps of model galaxies (seen edge on) 
exposed to a supersonic IGM wind, at T = 0.4 /i~-'^Gyr. The two 
panels show results for SPH (left) and VPH (right), respectively. 
Interestingly, the SPH galaxy is accelerated more strongly in the 
downstream direction, as a result of a larger effective ram pres- 
sure force due to slower stripping of its gaseous disk. The SPH 
simulation also shows signatures of the wind imprinted onto the 
morphology of the stellar disk, a behaviour that is not seen in 
VPH or AREPO. The two maps show the same spatial region of 
the wind tunnel and have an extension of 30 h~^kpc x 45 h~^]<ipc. 



galaxy, they eventually violate the distance condition and 
are then counted as stripped gas. 

In Figure [9] we show the ratio of the remaining ISM gas 
mass to the initial mass a function of time, comparing wind 
tunnel simulations carried out with the three different nu- 
merical techniques. Clearly, the mass loss in SPH is smallest 
overall, with about 10% of the initial mass remaining after 



1 Gyr, whereas at this time the ISM of the VPH and AREPO 
runs is almost stripped completely. Note that this measure 
of the mass loss counts most of the dense gas blobs that 
stay coherent in SPH as lost gas. If we were to measure the 
amount of "non-dispersed" gas instead, the difference would 
be considerably larger because the gas that is stripped in 
SPH from the disk is largely not dispersed, unlike in the 
other two approaches (see below). The relative similarity of 
the overall shape of the gas mass evolution in Figure [9] is 
therefore not simply a time delay in SPH relative to the 
other methods, even though at intermediate times the mass 
loss rates are similar. 

It is interesting to note that the stripped dense gas in 
the SPH simulation, and to a lesser extent in the VPH run, 
can continue to form stars at some level. This is particu- 
larly evident in Figure IIOI which shows maps of the star 
formation rate density in an edge-on projection, at time 
t = 0.4/i~^Gyr. This highlights that in the case of the 
SPH simulations the stripped gas clumps are not only very 
concentrated, they also continue to form stars. This phe- 
nomenon is still present in VPH, albeit at a much weaker 
level and with decreasing star formation rate over distance. 
In contrast, AREPO does not show such a behaviour. 

In SPH, there are at least two effects that can help to 
explain the enhanced star formation in the stripped clumps. 
First, it is known that spurious surface tension forces in SPH 
exist that will slightly compress the gas of the clumps and 
keep them coherent. Second, there is the issue of a potential 
artificial suppression of gas mixing in the turbulent wake 
behind the galaxy. Both SPH and VPH, by construction, 
prevent that the low entropy gas of the ICM is mixed at 
the particle level with the high entropy gas of the impacting 
wind. This imposed Lagrangian behaviour ignores any small- 
scale mixing processes that could happen on scales below 
the spatial resolution limit. In contrast, the mesh-based ap- 
proach of AREPO implicitly allows for such processes. Here, 
the mixing manifests itself through mass exchanges between 
the cells, giving rise to a much more diffuse wake behind the 
galaxy (left column in Figure [8]) where no gas remains that 
is still dense enough to support star formation. In the mesh- 
based approach, one may have in fact an opposite problem 
compared to the particle schemes. Here the need to advect 
the fluid over mesh interfaces can easily lead to excessive 
numerical mixing. However, the moving- mesh technique of 
AREPO should reduce such errors significantly compared to 
more traditional Eulerian methods. 

Another view on the different behaviour with respect 
to stripping is provided by Figure [TTl which shows the time 
evolution of the sum of the star forming gas mass (which 
is just the total gas mass above the star formation thresh- 
old) plus the stellar mass formed since the beginning of the 
simulation. Interestingly, this quantity is nearly constant in 
time in the case of SPH, showing that essentially none of the 
initially dense and star forming gas is dragged to low enough 
density to stop star formation. This is consistent with our 
earlier findings, but also makes it quantitatively evident that 
essentially no mixing of this low entropy gas with higher en- 
tropy material occurs. In contrast, the sum of star forming 
gas mass plus stellar mass declines significantly in AREPO, 
by more than a factor of 2 over the timescale of 1 /i~^Gyr 
that is shown here. The results for VPH are intermediate, 
showing some evidence for gas moving to higher entropy or 
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Figure 13. Resolution study of the mass loss of a galaxy in a wind tunnel. The left panel measures the gas remaining in the ISM (defined 
through Eq. [2]) as a function of time, for simulations carried out with VPH, SPH, and AREPO, as labelled. In each case, four resolutions 
are shown, corresponding to 4000 gas particles (very low resolution, dotted-dashed),to 20000 gas particles (low resolution, dotted), 80000 
gas particles (intermediate resolution, dashed), and 320000 particles (high resolution, solid) in the initial galaxy model. The panel on the 
right compares the sum of the remaining star forming gas plus the newly formed stellar mass in the whole simulation box as a function 
of time in VPH and SPH, using the same simulation set. Here a quite good convergence can be observed for the different particle-based 
techniques, but SPH is offset relative to VPH and has a significantly smaller amount of stripped gas. 



lower densities during the stripping process and thus reduc- 
ing the amount of total gas available for star formation. 

Interestingly, the discrepancies in the stripping rates 
also lead to different accelerations of the whole galaxies, 
because they experience unequal ram pressure forces as a 
result of different effective areas of the remaining gas disks. 
Since the gas in the disk is gravitationally coupled to the 
stellar disk and the dark matter halo, it is important to 
note that the wind not only strips the gas component of the 
galaxy, it also accelerates the whole galaxy. In Figure I12[ we 
show maps of the projected stellar densities of the stellar 
disks in SPH and VPH at an identical time, corresponding 
to T = 0.4/i~^Gyr. Clearly, the whole SPH galaxy has been 
pushed more into the downstream direction, demonstrating 
that it has experienced a higher effective ram pressure, as a 
consequence of the reduced stripping rate of its disk. In real 
galaxy clusters, this will in turn lead to a larger dynamical 
friction force and a somewhat faster decay of the orbit of the 
galaxy. An additional effect is that in the SPH case we see 
small features emanating in the downstream direction from 
the stellar disk. Figure [12] shows two small humps in the 
stellar disk, as well as a faint stellar blob at some distance 
further downstream. These are in fact stars of the original 
stellar disk that have been gravitationally pulled out of the 
disk by dense gas. Both effects are absent in the equivalent 
VPH and AREPO calculations, emphasising again how much 
more (artificial) "coherence" the dense gas phase exhibits in 
the case of SPH, where comparatively massive gas clumps 
are removed from the disk that stay coherent afterwards. 



Code resolutions soft, baryons softening DM 

VPH 4k, 20k, 80k, 320k 0.05 /j-^kpc 0.1 /z-^kpc 

SPH 4k, 20k, 80k, 320k 0.05 /i-^kpc 0.1 /i-^kpc 

AREPO 4k, 20k, 80k, 320k 0.05 /i-^kpc 0.1 /i-^kpc 



Table 2. Gravitational softening lengths for simulations in sec- 
tion0] The column "resolutions" denotes the initial number of gas 
particles in the disk where "k" denotes one thousand. The follow- 
ing columns show the gravitational softening lengths for baryons, 
that is for gas and stellar particles. The column on the right gives 
these values for the collisionless DM halo particles. 



4.3 Dependence on resolution and artificial 
viscosity 

In order to examine the quantitative robustness of our wind 
tunnel results with respect to numerical resolution, we re- 
peated our simulations with SPH, VPH and AREPO using 
both a lower and a higher resolution than in our default runs 
discussed thus far, where the galaxy is resolved with 80000 
gas particles in the initial gas disk. With respect to our de- 
fault resolution, our low resolution simulations have 4 times 
fewer particles in each component, whereas the high resolu- 
tion runs have 4 times mores particles, as denoted in Table (2] 
We here use gravitational softening lengths independent of 
resolution. 

In Figure [131 we show results of our resolution study 
for the gas stripping out of the ISM as a function of time, 
and for the sum of remaining star forming gas and newly 
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Figure 14. Dependence of the wind tunnel stripping efficiency in 
VPH (blue lines) and SPH (black lines) on the artificial viscosity. 
We show the loss of star-forming gas mass (defined as in Eq. ^ of 
our default galaxy model as a function of time, for high artificial 
viscosity a = 1 (solid lines), for intermediate strength o; = 0.5 
(dotted lines) and for low artificial viscosity a = 0.25 (dashed 
lines). 



formed stellar mass. We find that the general trends remain 
very similar at all three resolution levels, i.e., gas is stripped 
the fastest in AREPO and slowest in SPH, with VPH taking 
on an intermediate role. Interestingly, in the two particle 
codes, the stripping proceeds more slowly when the reso- 
lution is poorer, whereas in AREPO the opposite trend is 
observed, here the stripping tends to accelerate slightly for 
better resolution. We hence find that the results become 
more similar for better resolution. In particular, in the high 
resolution case, VPH is in fact quite close to AREPO. When 
the sum of the star forming gas mass and the newly formed 
stars is considered, we see that this quantity is very robust 
for SPH and VPH as a function of resolution, with the off- 
set probably simply reflecting different systematics related 
to the density estimation at the interface between gas disk 
and surrounding medium, as in Figure[7l On the other hand, 
the relative constancy of this quantity at late times reflects 
the inherent inability of these particle schemes to mix low 
entropy gas with gas of much higher entropy. 

Both SPH and VPH rely on an artificial viscosity to 
capture shock waves and to damp out small scale numerical 
noise. The choice of the artificial viscosity parametrisation 
and the associated strength of the viscous forces are crucial 
for the overall robustness and accuracy of the scheme. Using 
a large viscosity leads to more robust shock capturing, but 
on the other hand it may produce substantial viscous effects 
in region s where the gas really ought to flow without dis- 
sipation (|Cullen k Dehnenll201Ql ). In particular, a too high 
viscosity may suppress the growth of fluid instabilities that 
are important in gas stripping processes. We have therefore 
checked how our wind tunnel results for VPH depend on the 
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Figure 15. Dependence of the stripping efficiency in VPH on 
the strength of mesh shape correction forces, which consist of 
small, non-dissipative forces that tend to steer mesh cells towards 
a 'rounder' shape in case they have become highly asymmetric. 
We show the loss of star-forming gas mass (defined as in Eq. [2]) of 
our default galaxy model as a function of time, for our standard 
choice of /3o = 0.2, /5i = 0.01 (solid line) for parametrisation the 
shape correction forces. We also show results for two different set- 
tings, for /3o = 0.6, ^1 = 0.03 (dashed line) and for the somewhat 
extreme choice of /3o = 1-8, /3i = 0.09 (dotted line), respectively. 



viscosity setting. To this end we have reduced the viscosity 
parameter from our default value of a = 1 (high viscosity) 
to a = 0.5 (intermediate) and a = 0.25 (low artificial vis- 
cosity) . In Figure [M] we compare the corresponding results 
for the stripping as a function of time in our low-resolution 
wind tunnel set-up. We see that for reduced viscosity the gas 
is actually stripped a bit faster in VPH than for the high vis- 
cosity setting. This is consistent with our expectation that 
the numerical viscosity will tend to damp small-scale fluid 
instabilities and turbulence, which as a side effect reduces 
the stripping rate. 

Another numerical nuisance parameter in VPH lies in 
the strength of the shape correction forces that can be intro- 
duced into the technique in order to regularise the mesh and 
en courage a quasi- r egular particle distribution. As explained 
in lHeB Springel (|2010[ ), this is accomplished by adding a 
small term to the Lagrangian that penalises large aspect- 
ratio distortions of cells and large offsets between the gener- 
ating point of a cell and its centre-of-mass. While not strictly 
necessary for VPH to work, we have found the technique to 
be considerably less noisy in certain conditions when such 
weak shape correction forces are used. In our standard im- 
plementation, their strength is controlled by two dimension- 
less p arameters /3o — 0.2 and /3i = 0.01 (see iHeB k Sprii^ 
|2Q1Q| , for the exact definition of these parameters). In Fig- 
ure [151 we show the dependence of the stripping results 
for VPH when instead stronger shape correction forces de- 
scribed by f3o = 0.6, /3i = 0.03 or even /3o = 1.8, /3i = 0.09 
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Figure 16. Properties of the 'headwind' encountered by the in-faUing galaxy. The velocity difference is calculated between the centre- 
of-mass of the gas in the galaxy (defined according to Eqn. O and the environment, which is taken as a sector of a spherical shell with 
an opening angle of 30° in the direction of the velocity vector and at a distance of [45, 99] /i~-'^kpc with respect to the galaxy centre. All 
simulations have been conducted with cooling and star formation (the cited masses and times carry a factor of in their units). In 
each panel, we show results for simulations with VPH (blue), SPH (black), and AREPO (red), comparing from top left to bottom right 
the relative velocity, density, pressure and temperature of the wind encountered by the galaxy as a function of time. 
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are invoked. In both cases, this has only a very weak influ- 
ence on the results, as desired. 



5 STRIPPING OF A GALAXY DURING 
CLUSTER IN-FALL 

In order to complement our wind-tunnel experiments carried 
out in the previous section with more realistic setups, we 
here want to conduct a few simulations where galaxy mod- 
els are in-falling into live models of galaxy clusters. This ap- 
proach includes a full treatment of gravity as well as a mod- 
elling of cooling and star formation. In particular, it accounts 
correctly for the tidal effects of a galaxy travelling within the 
cluster potential. The latter inevitably changes the structure 
of the galaxy, most prominently by reducing its dark mat- 
ter mass through tidal truncation, which in turn changes 
the conditions under which the hydrodynamical processes 
occur. 



Besides looking at the stripping of the gas, we will 
also compare how the star formation rates in the in- 
falling galaxy decline in our different simulations, since 
we expect that different numeric al ram pressure will 
strongly affect this quantity (e.g. Kronberger et al.] l2008l : 
iKapferer eralll2009l ). We note that a number of recent stud- 
ies b oth with ^ rid-based codes (Roediger & Briiggen 2 0071: 
lapic hino k, Nie mever 2008) and SPH codes (Jachym et al.l 
I2OO7I : iMcCarthv et al.l l2008) have begun to investigate this 
question in detail. Our focus here will be much more limited 
though, and only be concerned with a comparison of our 
new VPH scheme relative to SPH and AREPO. 



5.1 Setup of galaxy-cluster interaction simulations 

For definiteness, we consider a parabolic encounter of 
a 10^^ /i~^M0 disk galaxy with a small M200 = 5 x 
lO^^/i~^M0 galaxy cluster. The galaxy is constructed as a 
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Figure 17. Visual comparison of the spatial distribution of stripped gas particles (large blue dots) from the ISM of an in-falling galaxy 
simulated with VPH (top row) and SPH (bottom row) at t = 0.3 Gyr (left column), t = 0.6 Gyr (middle column) and t = 0.9 Gyr (right 
column) (the cited masses and times carry a factor of in their units). The green curve traces out the orbit of the galaxy through the 
cluster up to the time of each individual panel, as labelled. The gas particles of the target cluster's IGM are shown as yellow dots. 



compound system as described in Section [3] but scaled down 
to fewer particles than in our earlier simulations. This was 
necessary to reduce the computational cost to an acceptable 
level, given that in this setup we need to represent the much 
larger cluster with the same mass resolution as the galaxy in 
order to avoid numerical problems in SPH (Ott & Schnette^ 
[2003l : lRead k Havfieldll201ir ). In our default set-up, we have 
therefore chosen 4000 gas particles for the ISM of the in- 
falling galaxy, and 10^ gas particles of identical mass for the 
IGM of the galaxy cluster. This low resolution will of course 
limit our ability to properly resolve hydrodynamical insta- 
bilities. We note however that similar and often still lower 
resolution is routinely used in cosmological simulations of 
galaxy formation, so our results are directly representative 
of such studies. They in any case allow us the extend our 
previous tests in Section |4] to lower resolution and investi- 
gate how reliably they can be assumed to carry over to more 
complicated situations. 

For the cluster, we adopt a gas-fraction of 17% and 
a NFW concentration of c = 2.0. A detailed list of main 
numerical parameters adopted for our simulations in this 
section is given in Table O The trajectory of the encounter 
starts at a distance of 700/i~^kpc, just outside of the virial 
radius, and follows a parabolic orbit where the minimal dis- 
tance of the two objects would be 50/i~^kpc if they were 
point masses. Relative to the orbital plane, the galaxy's 
disk is tilted by — 30° and th en turned by = 30° in 
the azimuth fsee lPuc et ahlEoOOl . for a sketch of the orbital 



geometry). We note that we can expect to be able to draw 
quite general conclusions from a single choice of inclination 
angles since the orbital geometry has b een shown to have no 
significant effect on the gas stripping (|Roediger &; Briiggenl 
l2006|). 



5.2 Properties of the head wind 

In Figure [ini we show the properties of the wind encountered 
by the in-falling galaxy as a function of time. We give the 
time evolution of the galaxy's environment as found in a 
sector with opening angle 30° of a spherical shell in the 
radial range 45 h~^kpc < r < 99/i~^kpc centred around the 
direction the galaxy is heading to. The velocity difference 
'i^headwind showu in the fist panel of Figure [16] is computed 
with respect to the centre of mass of the galaxy. The figure 
shows that density, pressure and temperature of the head 
wind reach their peak at the pericentre, as expected. 

We include results for VPH, SPH and AREPO in Fig.fTBl 
and as the comparison shows, there is generally very good 
agreement between the three different simulation techniques. 
All three runs show a strong rise of density and pressure as 
the galaxy approaches and passes pericentre, while the tem- 
perature is roughly constant, reflecting the nearly isother- 
mal conditions in the cluster gas. There are no significant 
deviations in the orbit of the galaxy between the different 
methods, hence any difference in the evolution of the galax- 
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Code 


softening baryons 


softening DM 


TsN 


Ao 


regularisation ^q/i 


VPH 


0.25 /i-^kpc 


0.5 h-^kpc 


3 X 10^ 


3000 


/3o = 0.6, /3i = 0.03 


SPH 


0.25 h-^kpc 


0.5 h-^kpc 


3 X 10^ 


3000 




AREPO 


0.25 h-^kpc 


0.5 h-^kpc 


3 X 10^ 


3000 





Table 3. Gravitational softening lengths for simulations in section [5] The columns show the gravitational softening lengths for baryons, 
that is for gas and stellar particles. The next column shows these values for collisionless DM halo particles. The following columns specify 
the parameters of the supernova feedback model, where TgN is the effective temperature and Aq the evaporation efficiency (using the 
notation of Springel Sz Hernquist 2003). The last column gives the parameters used for the Voronoi cell shape correction scheme in VPH 
(in the notation of .Hefi &: SDringeL.201ft) . 




o 
O 



0.0 0.2 



0.4 0.6 
t [Gyr/h] 



0.8 



SPH 
VPH 
Arepo 




0.0 0.2 



0.4 0.6 

[ h-^ Gyr ] 



0.8 



Figure 18. Gas mass as a function of time still bound to a galaxy 
that is in-falling into a galaxy cluster. We include simulation re- 
sults for SPH (black), VPH (red) and AREPO (blue). 



Figure 19. The stellar mass formed within a distance of 
15 h~^kpc from the galaxy centre as a function of time. We in- 
clude simulation results for SPH (black), VPH (blue) and AREPO 
(red). 



ies can only arise from differences in the hydrodynamical 
treatment of the interaction of galaxy and cluster gas. 

5.3 Gas stripping and star formation truncation 

In order to define whether a gas particle or cell still belongs 
to the galaxy we use the condition: 

p > 4 X 10"^^ g cm"^ and T < 10^ K. (3) 

Furthermore, we additionally require that the separation r 
of particles from the centre-of-mass of the remaining grav- 
itationally bound dark matter of the galaxy is less than 
20 h~^kpc. We note that with this definition the galaxy can 
in principle also accrete new gas from the cluster which may 
then cool onto the ISM and help to provide fuel for star 
formation. 

In Figure \17\ we show a visual comparison of the gas 
stripping in VPH and SPH simulations, where it is readily 
apparent that the gas mass loss proceeds much slower in 
SPH than in VPH. This is confirmed by quantitative mea- 
surements shown in Figure [TSl which include results for gas 
stripping in three different simulations of the encounter of a 



galaxy with the cluster ICM. Interestingly, VPH looses gas 
here even somewhat faster than AREPO, but both methods 
yield a substantially faster loss of gas than SPH. While at 
time ~ 1 Gyr the VPH simulation has lost all of the gas, the 
galaxy in SPH still retains half of it. 

This substantial difference is corroborated by the be- 
haviour of the star formation rate in the galaxies, shown in 
cumulative form in Figure [19] Here SPH shows the slowest 
decline overall, while both VPH and AREPO lead to a rapid 
termination of star formation, which implies a quick red- 
dening of the galaxies. It is a well known problem in galaxy 
formation to understand the colours of cluster galaxies in 
detail. Semi-analytic models of galaxy formation have com- 
monly predicted a rather quick truncation of star forma- 
tion upon cluster in-fall, yielding cluster populations that 
are actually too red (e.g. 'Weinm ann et al.ll20Q6l V It can be 
expected, due to the faster stripping, that our VPH and 
AREPO simulations suffer from the same problem, poten- 
tially making it even more severe. This may indicate that 
the stripping efficiencies in all of the simulations are sub- 
stantially too large, probably because the ISM is still under- 
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Figure 20. Gas fraction as a function of distance to the cluster 
centre for two representative (sub)halos that fall into the cluster. 
We show results for a substructure with 4.5 k DM-particles (top 
panel) and 1.8 k DM-particles (bottom) panel, comparing results 
for SPH (black), VPH (blue) and AREPO (red). We find a clear 
difference in the final phase of the stripping process. In the SPH 
simulation, the substructure keeps a higher gas content until it is 
disrupted. 



resolved and appears as a homogeneous dense phase instead 
of being resolved into a true multi-phase medium. The latter 
would make the medium more porous, allowing very dense 
clouds of gas to resist the ICM wind for a longer time and 
to stay in the in-falling galaxy. 



6 COSMOLOGICAL CLUSTER SIMULATIONS 

We now turn to results for fully cosmological simulations of 
cluster formation. We first carry out a comparison of the 
gas content of satellite systems in 'zoom' simulations of the 
formation of rich galaxy clusters, carried out with VPH, 
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Figure 21. Gas fraction of substructures as a function of their 
distance to the cluster centre. We compare results for SPH 
(black with green shade), VPH (blue) and AREPO (red with 
orange shadow), including only subhalos more massive than 
M > 3 X IQ-'^-'^ h~^M(7) to exclude poorly resolved small structures. 
The shaded areas indicate one standard deviation from counting 
statistics. 



SPH and AREPO. We here primarily want to check whether 
the differences we have observed in the stripping of dense 
ISM gas out of galaxies in our earlier more idealised simula- 
tions manifest themselves also in non-radiative simulations 
without star formation, where the density contrast is much 
smaller. Secondly, we study the well known Santa Barbara 
cluster o f iFrenk et al.l (|l999l l in order to investigate whether 
VPH produces a higher entropy in the cluster centre com- 
pared to SPH, which would then make it closer to the results 
of mesh codes that have been applied to this problem. We 
also use this cluster in order to assess the numerical conver- 
gence of VPH for the properties of the intra-cluster gas. 



6.1 Gas stripping in non- radiative zoom 
simulations of galaxy clusters 

In order to simulate the formation of rich galaxy clusters 
in the ACDM cosmology, w e extract a massive halo from 
the Millennium Simulation dSp ringel et"all l2QQ5l ). and re- 
simulate this cluster with the addition of gas. To this end 
we trace the particles that make up the cluster back to the 
original initial conditions at 2; = 127, thereby finding the 
Lagrangian region out of which the cluster has formed. This 
region is then populated both with dark matter and gas par- 
ticles, which are perturbed with the original displacement 
field. We can also increase the resolution compared to that 
in the original simulation if desired, in which case additional 
small-scale fluctuation power is added in the region between 
the old and new Nyquist frequencies. Further away from 
the region that holds the cluster material and its immediate 
surroundings, the resolution is reduced by combining parti- 
cles into progressively heavier 'boundary' particles. In this 
way, the resolution gradually decUnes with distance from the 
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cluster while the gravitational tidal field that influences its 
formation is still accurately determined. With this standard 
'zoom' technique, the computational cost is concentrated in 
the small region of interest, allowing high resolution simula- 
tions of individual objects in comparably short time. 

For definiteness, we have picked a cluster of virial mass 
M200 = 1.8 X lO^^/i"^M0, which we re-simulate with a 
baryon fraction = 0.045. The other cosmological param- 
eters are the same as in in the original Millennium run and 
are given by Qm = 0.25, Qa = 0.75, ag = 0.9, h = 0.73, and 
n = 1. We have initialised the re-simulations at 2; = 127 and 
evolved them with VPH, SPH and AREPO to the present 
epoch. We employed an identical gravitational softening 
length of 0.005 /i~^kpc for all particle types and hydrody- 
namical schemes, and treated the gas as a non-radiative mix 
of hydrogen and helium. 

We identify halos in the simulations by applying the 
FOF algorithm to the high-resolution dark matter particles. 
Each gas particle is assigned to the group in which its closest 
dark matter particle lies. W e then apply the SUBFIND algo- 
rithm (Springel et"alll200l ') to the particle groups we found 
in this way in order to decompose them into gravitationally 
bound (sub)groups. Using the IDs attached to each parti- 
cle we can track individual halos/subhalos as a function of 
time, and, in particular, study how the gas content of subha- 
los declines as a halo falls into the cluster. In order to reduce 
numerical noise in our simulation comparison, we however 
restrict our substructure selection to a sample where more 
than 60% of the DM particles can be found in every studied 
simulation run. 

As an example, we show in Figure [20] the evolution of 
the gas fraction of two different substructures as a func- 
tion of their distance to the cluster centre (the correspond- 
ing times are indicated by the redshift labels on the upper 
axis), comparing results for VPH, SPH and AREPO simula- 
tions. For both substructures, we find good agreement in the 
early in-fall phase, for distances larger than about ^ 3i?vir- 
At smaller separations, the influence of the cluster is how- 
ever very noticeable already, and the substructures loose gas 
quickly. In both of these examples, we identify a consider- 
ably larger gas loss in the VPH run than in SPH at the last 
time when we can still find the substructures before they 
are disrupted. This is consistent with our earlier findings for 
the stripping of the dense ISM gas. 

In Figure [21] we now consider all of the substructures 
around the cluster at the final time, and simply compare the 
gas mass fraction in substructures as a function of cluster- 
centric radius. Even though the results are a bit noisy, we 
find a clear statistical trend of a smaller gas mass fraction 
in the VPH run relative to SPH, especially in the region 
^ 1 — 3 i?vir, where the gas fraction declines rapidly. Interest- 
ingly, a simulation with AREPO for the same cluster initial 
conditions shows a bound gas fraction that is still smaller 
in-falling dark matter halos. We hence conclude that even 
at the level of non-radiative simulations there are already 
significant differences in the stripping of gas out of in-falling 
substructure, which also implies that the mixing of gas in 
massive halos can be expected to differ in the three exam- 
ined numerical techniques. 



6.2 The Santa Barbara Cluster 

The ^ Santa Barbara cluster comparison project' of 
iFrenk et al. 1 {1999) analysed the results of a large number 
of different cosmological hydrodynamical codes for the for- 
mation of a rich cluster of galaxies with non-radiative gas. 
The comparison involved both SPH codes and hydrodynam- 
ical mesh codes, and focused, in particular, on the resulting 
thermodynamic properties of the intra-cluster gas. An im- 
portant result that emerged from the study was that the 
different methods systematically disagree in the amount of 
entropy predicted for the central regions of the cluster, with 
the mesh-based approaches yielding consistently higher cen- 
tral entropy and correspondingly lower density than the SPH 
codes. 

The Santa Barbara (SB) cluster has become a standard 
test problem for cosmological hydrodynamic codes, with re- 
sults reported in numerous studies. Recently, some studies 
have suggested that the difference seen between the var- 
ious techniques is prima rily associated with differences i n 
the treatment of mixing ('Mi tchell et~alll2009l : IVazzall201ll ). 
which is suppressed in SPH by construction, and this may 
artificially lower the central entropy. Specifically, it has been 
suggested that this problem may be related t o a suppression 
of Raleigh- Taylor fluid instability in SPH (|Wadslev et al.l 
2008). If the difference is caused by the lack of mixing at 
the particle level (Tasker et al. 2008; Wadslev et al. 2008l l 
the VPH results should in principle agree with SPH. 

In the spirit of the original project of iFrenk et al.l 

([l99i), we here re-run the SB cluster with our new hydro- 
dynamical VPH method. We are especially interested in the 
question whether VPH differs in its predictions for the cen- 
tral cluster region compared with SPH, which perhaps could 
arise due to the different stripping efficiency of this tech- 
nique. We use the same initial conditions that have been 
used in the original SB cluster comparison project, where 
an Einstein-de-Sitter cosmological model with mean den- 
sity Qm = 1, Hubble constant Ho = 50 km s~^Mpc~^, and 
baryon fraction of = 0.1 was used. The initial conditions 
were constructed as a constrained realisation where a rich 
cluster corresponding to a 3a peak was imposed to form in 
the centre of a cubic box with a side-length of 64 Mpc, in an 
otherwise random realisation. 

We have simulated the SB cluster at a variety of differ- 
ent resolutions using VPH, ranging from 2 X 32^ to 2 X 256^ 
particles with parameters as denoted in Table [4] in order to 
see how well the method converges for the primary cluster 
properties. In Figure [22] we show the results of this conver- 
gence test, in terms of radial profiles for gas density, dark 
matter density, enclosed gas fraction, temperature, entropy, 
and specific X-ray emissivity. We find that the convergence 
is in general good for the outer profiles and the dark matter 
properties. Only in the very centre some interesting differ- 
ences can be observed. For the most part they can be simply 
understood as an imprint of the varying spatial and mass 
resolution across the sequence, where shortly before the res- 
olution limit is reached the lower resolution runs peel away 
from the converged result seen in the higher resolution sim- 
ulations. This is for example the case for the entropy and 
gas density profiles. 

It is clear however that VPH does not predict the forma- 
tion of a large constant entropy core region, unlike what is 
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Figure 22. Resolution study with VPH for the cluster of the cosmological "Santa Barbara Cluster Comparison project". Radial profiles 
of (from top left to bottom right) give gas-density, dark matter density, cumulative gas-fraction, gas-temperature, specific entropy and 
X-ray luminosity. All properties have been obtained as mass-weighted arithmetic averages within logarithmically spaced bins (the cited 
masses and times carry a factor of in their units). The different colours give VPH results with initial gas particle numbers of 32^ 
(black), 64^ (red), 128^ (blue) and 256^ (green), respectively. 
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^, ^1 


= 0.04 


VPH 
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"-•^kpc 
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^0 


= Oi 


^, ^1 


= 0.04 
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^, /3i 
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VPH 


2 X 256^ 


3 h~^kpc 
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20 K 


^0 


= Oi 


^, ^1 


= 0.04 



Table 4. List of parameters for simulations in section [6.21 The column "resolutions" denotes the initial particle numbers. The following 
columns show the gravitational softening lengths for gas and collisionless DM halo particles. The columns on the right specify parameters 
deviating from our standard choice, such as minimal gas temperature and slightly more aggressive values for VPH's mesh regularisation 
scheme. 



typically seen in mesh codes and consistent with the falling 
entropy profiles that are typically observed with SPH all 
the way to the centre of halos. In fact, in general our results 
obtained with VPH for the Santa Barbara cluster are in 
^oo d ag;reement with SPH results reported i n the literatur e 
(e.g. iFrenk et al.]ll999l : lAscasibar et al.ll2003l : ISpringelllioOsh , 
suggesting that it makes comparatively little difference for 
the bulk thermodynamic structures which particle-based 
method is used in non-radiative cosmological simulations. 
In particular, the different stripping efficiency alone and the 
slightly higher spatial resolving power of the Voronoi-based 
density estimate do not yield an entropy core as generally 
found in hydrodynamical mesh-codes. 



7 DISCUSSION AND CONCLUSIONS 

In this study, we have carried out a systematic comparison 
of the properties of our new Voronoi part icle hydrodynam- 
ics (VPH) method (jHefi Springell I2OIOI I with respect to 
standard SPH and moving-mesh hydrodynamics as imple- 
mented in the AREPO code. We have focused on stripping 
processes in galaxies and halos upon in-fall into galaxy clus- 
ters. Here, it is expected that the outer parts of gaseous 
disks are qu ickly removed due to ram pressure stripping 
(|Gunn &: Go tt 1972), but the subsequent more gradual gas 
loss sensitively depends on the ability of hydrodynamical 
codes to capture fluid instabilities occurring in shear flows 
around the galaxies. The recent flndings that SPH appears 
to exhibit severe inaccuracies in this regime has prompted us 
to develop the alternative VPH method. In the present pa- 
per, we have studied how this new technique compares with 
traditional SPH and the new moving-mesh code technique 
in a number of basic problems relevant for galaxy formation. 

To this end, we flrst compared results for isolated com- 
pound galaxy models, both in isolation and in wind tun- 
nels where they were exposed to a supersonic head wind. 
This set-up allowed relatively high-resolution simulations of 
wind-ISM interactions. Our simulations have revealed non- 
negligible differences in the rate at which dense ISM gas is 
stripped and dispersed, and in the appearance of this gas in 
the downstream part of the flow. SPH showed a lower strip- 
ping rate than both VPH and the mesh-code AREPO. As a 
result, the SPH galaxy also experienced the largest displace- 
ment due to the ram pressure of the impinging wind. Also, 
we were able to show that essentially none of the ISM gas in 
SPH and VPH could ever be transferred to much lower den- 
sity. Instead, if gas was stripped, it stayed in coherent dense 
blobs, where even star formation could continue. Further- 
more due to elevated pressure caused by the surface tension 



effect in SPH, the star formation in these stripped blobs re- 
mained at an unphysically high level. AREPO, in contrast, 
showed a rapid loss of gas out of the disk, which was further- 
more efficiently mixed with other gas, so that lower densities 
were reached quickly by the stripped gas and star formation 
was stopped. 

We followed up these simulations with numerical exper- 
iments where we dropped galaxy models in live cluster mod- 
els. Even though here the resolution was substantially lower, 
we obtained results in good qualitative agreement with our 
wind tunnel runs. Likewise, in non-radiative cosmological 
simulations of galaxy cluster formation, we followed indi- 
vidual subhalos as they fell into the forming cluster, find- 
ing again that the gas content of satellite systems declined 
slowest in SPH, while the stripping in VPH and especially 
in AREPO proceeded noticeably faster. 

Finally, we considered simulations of the Santa Barbara 
cluster, which has become an important test problem for 
evaluating cosmological hydrodynamical codes. While the 
VPH runs revealed a slightly elevated entropy compared 
to SPH at the smallest radii, they in general agreed quite 
well with SPH and in particular did not provide evidence at 
the highest resolution that an entropy core similar to those 
found in mesh-codes such as AREPO is formed. This is per- 
haps to be expected if the e ntropy core indeed primarily 
arises from mixing processes (jMitchefl et al.l 120091 ) that are 
largely absent in SPH and VPH, by construction. 

Overall, it thus appears that VPH offers some improve- 
ments over SPH without however changing its fundamental 
character. We argue that the most important origin of these 
differences lies in an improved gradient estimate in VPH 
compared to SPH. VPH is second-order accurate in the gra- 
dient estimates, i.e. a linear gradient is always reproduced 
exactl y, independent of the particle distribution (|SDringell 
l2010al i. In contrast, SPH has a so-called zero-th order error 
in its gradient estimate (e.g. Read et al. 2010). This in par- 
ticular means that even for equal pressure s for all par ticles 
the pressure force not necessarily vanishes ([AbellbOllh , and 
furthermore, the absolute size of the gradient error grows 
with the pressure itself. The gradient errors in SPH have 
also been linked to the creation of small-scale velocity noi se 
in studies of subsonic turbulence (jBauer k Springelll20Tlh . 

A good test problem for appreciating this difference in 
the gradi ent accuracy is the flow of a stable vortex as sug- 
gested bv lGresho k ChanI (|l990l ). In this 'triangular vortex 
problem' a fluid is set up with an azimuthal velocity proflle 
(see Appendix [Al for details), such that the pressure gra- 
dients counter the centrifugal force, and the vortex evolves 
in a time-independent, stable fashion. Figure [23] shows the 
radial velocity profile after the vortex has been evolved for 
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Figure 23. Radial profiles of the azimuthal velocity in evolved 
simulations (t = 1.0) of the Gresho vortex test, calculated at 
resolution 80x80 with different numerical techniques (SPH, VPH, 
and AREPO), as labelled. The small black dotes show individual 
velocities of simulation particles/cells, while the red dots give the 
averaged solution binned in radial annuli. The blue thick lines 
give the analytic solution that is realised in the initial conditions. 
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Figure 24. Convergence rate of three different numerical tech- 
niques for the Gresho vortex test. The black points refer to mea- 
surements of the LI error norm at t = 2.0 for the SPH method, 
with resolutions ranging from 20 x 20 to 640 x 640. The LI er- 
ror has been calculated on the basis of all single particles/cells. 
The blue dots give the corresponding measurements for our VPH 
method, while the red points are for the AREPO moving-mesh 
code. The dashed lines indicate power-laws with slopes —0.3, —0.7 
and -1.4 for SPH, VPH and AREPO, respectively. 



a time t = 1 with the VPH, SPH and AREPO codes in 2D, 
using a 80 x 80 Cartesian grid for the initial conditions in the 
domain [-0.5, 0.5]^ It can be seen clearly that SPH shows a 
much larger velocity scatter than the other two codes, and 
its solution has already degraded quite noticeably relative 
to VPH and AREPO. Especially in the inner domain, where 
the fluid rotates like a solid body, SPH deviates systemat- 
ically from the analytic solution. We note that part of this 
degr adation can be influenced by the artific ial viscosity set- 
ting ("Spri ngelll2010bl : lRead HavfieldllioTT ') , but if a higher 
viscosity is used to suppress the velocity noise it typically 
also leads to a faster viscous erosion of the vortex. In VPH, 
some velocity scatter is seen as well, but it is appreciably 
smaller than in SPH, which we interpret as a consequence 
of the more accurate gradient estimates in VPH. 

Arguably one of the best ways to quantify the accuracy 
of the different numerical techniques for this analytic test 
problem is to look at an objective error measure as a func- 
tion of resolution. In Figure I24[ we compare the Ll-error 
norm for the azimuthal streaming velocity of the Gresho 
test as a function of resolution for all three techniques. It 
is evident that ordinary SPH converges only very slowly, 
whereas VPH shows a considerably improved convergence 
rate. This directly demonstrates an important improvement 
brought about by VPH, which we can directly trace to better 
gradient estimates. The latter appears also as the primary 
reason for the better results for stripping we obtained with 
VPH compared to SPH, which is due to the more accurate 
treatment of contact discontinuities and the avoidance of the 
'gap' phenomenon of SPH. 
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Note that the difference in convergence rate also means 
that VPH will outperform SPH in terms of computational 
cost if very high accuracy needs to be achieved, provided 
its computational cost differs only by a constant factor of 
order unity which is indeed the case in our present imple- 
mentation. However, the Voronoi mesh construction adds 
substantial computational cost compared to ordinary "stan- 
dard" SPH, making the VPH code about a factor 3-4 slower 
for pure hydrodynamics when the same number of particles 
is used. This is mitigated to some extent if self-gravity is 
included (which is typically about as expensive or slightly 
more expensive than SPH-based hydrodynamics), reducing 
the difference to less than a factor of 2. We note that some 
alternative su ggestions to impro ve standard SPH, such as 
the scheme bv lRead et al.l (|2010l l which involves a dramatic 
increase of the number of smoothing neighbours, also re- 
quire an increased computational cost per resolution ele- 
ment. Which of these schemes is ultimately the most efficient 
one (in the sense of reaching a given accuracy goal with the 
smallest computational cost) is difficult to answer in general, 
and is in any case a problem-dependent and implementation- 
dependent question. 

According to Fig. El VPH still falls short of the bet- 
ter convergence rate obtained with the moving-mesh code 
AREPO (and similarly with fixed grid m esh codes such as 
ATHENA, see lStone^ et al.ll2008l : ISDringelll2Q10al ). The above 
discussion suggests that higher order density estimates com- 
bined with at least equally accurate gradient estimates are 
needed to improve on SPH and VPH in this respect. Some 
suggestions in this direction have recently been made (e.g. 
Read k Havfieldl l2Qlll : iMaron et al.1 l2Qlll : iMcNallv et al.1 
201lh . and it will be interesting to see whether they can 
successfully yield significant accuracy improvements in cos- 
mological applications such as those discussed here. 
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APPENDIX A: GRESHO VORTEX TEST 

In the test problem of a stable vortex as suggested by 
iGresho k ChaJ (Il990l ) the azimuthal velocity profile has the 
form 
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Figure Al. Particle distribution and resulting density measure- 
ments in the case of perfectly circular particle motions. The four 
panels on top illustrate the distortion and shearing of the ini- 
tial Cartesian grid if the points would move on perfectly circu- 
lar paths with their initial azimuthal velocity. Such orbits would 
be obtained if the numerical pressure gradient would always 
agree with the analytically expected gradient. The bottom panel 
gives the density estimates obtained with VPH (blue points) and 
SPH (black circles) for the particle distribution obtained at time 
t = 0.8 (middle right panel). 



in a gas of constant density equal to p = 1 and adiabatic 
index of 7 = 5/3. The pressure is chosen as 

5 + 25/2r2 for ^ r < 0.2 
9 + 25/2r2- 

20r + 41n(r/0.2) for 0.2 ^ r < 0.4 

3 + 41n2 for r ^ 0.4 



P(r) = <^ 



(A2) 



so that the pressure gradients balance the centrifugal force. 

In Figure [23l we showed the radial velocity profile after 
the vortex has been evolved for a time t — 1 with the VPH, 
SPH and AREPO codes in 2D, using a 80 x 80 Cartesian grid 
for the initial conditions in the domain [—0.5, 0.5]^. Interest- 
ingly, one can also clearly see that the velocity noise in VPH 
is still lower in the solid body rotation part of the vortex, 
for r < 0.2. Here the initial pressure gradients are correct 
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in VPH, and they will in principle stay that way during the 
motion of the solid body part. However, the strongly shear- 
ing part between 0.2 < r < 0.4 will eventually disturb the 
evolution - and this can in some sense be blamed on the 
density estimate. 

To see this, let us imagine for the moment that the fluid 
particles would all move under the exact analytic pressure 
gradient, then they would move on circular orbits with con- 
stant angular velocity. The resulting time evolution of the 
point set is shown in the top four panels of Figure [ATI to- 
gether with the resulting density estimates at time t = 0.8 
for VPH and SPH in the lower panel. The shearing motion 
creates substantial irregularities in the particle distribution, 
manifesting itself in significant fluctuations of the density 
estimates around the desired background density of p = 1. 
These fluctuations tend to be roughly of comparable mag- 
nitude for the kernel-based density estimates in SPH and 
the Voronoi-based density estimates in VPH. In any case, 
the density fluctuations of course cause pressure fluctuation 
that necessarily prevent that the particles move on exactly 
circular orbits. Instead, additional motions are created that 
locally restore the pressure equilibrium. Motions that do not 
agree with the analytic characteristics are here needed to 
mitigate pressure fluctuations and to prevent that density 
errors as large as implied by the analytic trajectories occur 
in the first place. In contrast to SPH and VPH, in AREPO, 
the density values can stay correct because this code calcu- 
lates mass exchanges between cells in such a way that the 
density of cells stays (nearly exactly) constant, despite the 
strong shearing of the tessellation. Combined with the ac- 
curate gradient operator in this code (which is the same as 
in VPH), this avoids spurious contributions to the pressure 
gradients, leading to higher accuracy in the final solution. 

A further layer of complexity is added when the ambigu- 
ous role of the shape correction forces in VPH is considered. 
On one hand they induce deviations of particle motions from 
circular orbits in order to make cells more regular, which can 
degrade the accuracy with which the analytic solution is fol- 
lowed. On the other hand they maintain regularity of the cell 
configuration, which improves the robustness of the scheme 
and helps to ensure ac curate density and gradient estimates 
dHefi & SDringel'l2010|). It is well possible that an improved 
cell regularisation scheme is able to improve VPH's perfor- 
mance for the Gresho problem, something we leave for future 
investigations. 
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